1 Sinusoidal Model

Goal: to fit a simple sinusoidal^[1] model to the sunspots data.

Let's consider the expression of sinusoidal model: $$s(t)=A\cos (2\pi ft)+B\sin (2\pi ft),A=R\cos \varphi ,B=-R\sin \varphi .$$
For us, usually ${\displaystyle t=1,2,\cdots ,n}$. In this case, we can restrict the frequency ${\displaystyle f\in [0,\frac{1}{2}]}$.^[2]

Based on the remark,

When ${\displaystyle f=0}$, ${\displaystyle s(t)}$ is constant; when ${\displaystyle f=\frac{1}{2}}$, ${\displaystyle s(t)=R\cos (\pi t+\varphi )=(-1{)}^{t}R\cos (\varphi )}$.

2 Frequentist Inference

Like in linear models, frequentist calculate the MLE: $$\begin{array}{rl}& \prod _{i=1}^n\frac{1}{\sqrt{2\pi }\sigma }\exp [-\frac{(y_t-{\beta }_0-{\beta }_1\cos 2\pi ft-{\beta }_2\sin 2\pi ft{)}^2}{2{\sigma }^2}]\\ \propto & {\sigma }^{-n}\exp [-\frac{1}{2{\sigma }^2}\sum _{t=1}^n(y_t-{\beta }_0-{\beta }_1\cos 2\pi ft-{\beta }_2\sin 2\pi ft{)}^2]\\ =& {\sigma }^{-n}\exp [-\frac{1}{2{\sigma }^2}S({\beta }_0,{\beta }_1,{\beta }_2,f)].\end{array}$$
So the problem is to find ${\displaystyle ({\hat{\beta }}_0,{\hat{\beta }}_1,{\hat{\beta }}_2,\hat{f})}$ that minimizes $$\begin{array}{rl}S({\beta }_0,{\beta }_1,{\beta }_2,f)& =\sum _{t=1}^n(y_t-{\beta }_0-{\beta }_1\cos 2\pi ft-{\beta }_2\sin 2\pi ft{)}^2\\ \text{(2.1)}& & =||y-X_f\beta |{|}^2.\end{array}$$ We first fix ${\displaystyle f}$, so this is just a linear regression. We have $$\hat{\beta }(f)=(X_f^{\mathrm{T}}{X}_f{)}^{-1}{X}_f^{\mathrm{T}}y,X_f=\left[\begin{array}{ccc}1& \cos 2\pi f(1)& \sin 2\pi f(1)\\ ⋮& ⋮& ⋮\\ 1& \cos 2\pi f(n)& \sin 2\pi f(n)\end{array}\right].$$ Then $$\hat{f}=\arg \underset{f}{min}S(\hat{\beta }(f),f),\beta =({\beta }_0,{\beta }_1,{\beta }_2{)}^{\mathrm{T}}.$$

3 Bayesian Inference

The posterior is $${\sigma }^{-n}\exp [-\frac{S(\beta ,f)}{2{\sigma }^2}]\mathbf{1}\{0\le f\le \frac{1}{2}\},$$ and assume prior ${\displaystyle {\beta }_0,{\beta }_1,{\beta }_2,\log \sigma \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}}{\sim }\mathrm{Uniform}[-C,C]}$, ${\displaystyle f\sim \mathrm{Uniform}[0,\frac{1}{2}]}$.
Then^[3] $$f_{\beta ,\sigma ,f|data}={\sigma }^{-n-1}\exp [-\frac{S(\beta ,f)}{2{\sigma }^2}]\mathbf{1}\{0\le f\le \frac{1}{2}\}\mathbf{1}\{-C<{\beta }_0,{\beta }_1,{\beta }_2,\log \sigma <C\},$$ then the posterior density of ${\displaystyle f}$ $$\begin{array}{r}\propto \iint {\sigma }^{-n-1}\exp [-\frac{S(\beta ,f)}{2{\sigma }^2}]\mathrm{d}\beta \mathrm{d}\sigma .\end{array}$$

By Pythagorean identity, $$S(\beta ,f)=||y-X_f\beta |{|}^2=S({\hat{\beta }}_f,f)+(\beta -{\hat{\beta }}_f{)}^{\mathrm{T}}{X}_f^{\mathrm{T}}{X}_f(\beta -{\hat{\beta }}_f).$$
So further on $$\begin{array}{rl}\propto & \iint {\sigma }^{-n-1}\exp (-\frac{S({\hat{\beta }}_f,f)}{2{\sigma }^2})[-\frac{(\beta -{\hat{\beta }}_f{)}^{\mathrm{T}}{X}_f^{\mathrm{T}}{X}_f(\beta -{\hat{\beta }}_f)}{2{\sigma }^2}]\mathrm{d}\beta \mathrm{d}\sigma \\ \\ =& \int {\sigma }^{-n-1}\exp (-\frac{S({\hat{\beta }}_f,f)}{2{\sigma }^2})(2\pi {)}^{\frac{p}{2}}\sqrt{det({\sigma }^2(X_f^{\mathrm{T}}{X}_f{)}^{-1})}\\ \propto & \int {\sigma }^{-n-1}\exp (-\frac{S({\hat{\beta }}_f,f)}{2{\sigma }^2}){\sigma }^p|X_f^{\mathrm{T}}{X}_f{|}^{-\frac{1}{2}}\mathrm{d}\sigma \\ =& |X_f^{\mathrm{T}}{X}_f{|}^{-\frac{1}{2}}\int {\sigma }^{-n+p-1}\exp (-\frac{S({\hat{\beta }}_f,f)}{2{\sigma }^2})\mathrm{d}\sigma \\ =& {\int }_0^{\mathrm{\infty }}{t}^{-n+p-1}(S({\hat{\beta }}_f,f){)}^{-\frac{n-p}{2}}\exp (-\frac{1}{2t^2})\mathrm{d}t\\ \propto & (S({\hat{\beta }}_f,f){)}^{-\frac{n-p}{2}}|X_f^{\mathrm{T}}{X}_f{|}^{-\frac{1}{2}}.\end{array}$$
Compare with linear regression: $$\mathrm{Posterior}\propto (S(\beta ){)}^{-\frac{n}{2}}.$$

4 Fourier Frequency

${\displaystyle \mathrm{RSS}(f)}$ (${\displaystyle S}$ above measures how well the sinusoid at frequency ${\displaystyle f}$ fits the data). So we should take a grid of values for ${\displaystyle f}$, and then compute ${\displaystyle \mathrm{RSS}(f)}$ for each grid point.

Common Choice of Grid for ${\displaystyle f}$: (notate as ${\displaystyle G}$)

• ${\displaystyle n}$ is even, ${\displaystyle 0,\frac{1}{n},\frac{2}{n},\cdots ,\frac{1}{2}}$,
• ${\displaystyle n}$ is odd, ${\displaystyle 0,\frac{1}{n},\frac{2}{n},\cdots ,\frac{n-1}{2n}}$.

Note that ${\displaystyle n}$ is the size of data.

So our common choice inside ${\displaystyle G}$ are all Fourier frequencies lying in ${\displaystyle [0,\frac{1}{2}]}$.

When ${\displaystyle f\in G}$, plug in ${\displaystyle {\hat{\beta }}_f=(X_f^{\mathrm{T}}{X}_f{)}^{-1}{X}_t^{\mathrm{T}}y}$, $$\begin{array}{rl}\mathrm{RSS}(f)=& ||y-X_f{\hat{\beta }}_f|{|}^2=(y-X_f{\hat{\beta }}_f{)}^{\mathrm{T}}(y-X_f{\hat{\beta }}_f)\\ =& y^{\mathrm{T}}y-{\hat{\beta }}_f^{\mathrm{T}}{X}_f^{\mathrm{T}}y-y^{\mathrm{T}}{X}_f{\hat{\beta }}_f+{\hat{\beta }}_f^{\mathrm{T}}{X}_f^{\mathrm{T}}{X}_f{\hat{\beta }}_f\\ =& y^{\mathrm{T}}y-y^{\mathrm{T}}(X_f^{\mathrm{T}}{X}_f{)}^{-1}{X}_f^{\mathrm{T}}y-y^{\mathrm{T}}{X}_f(X_f^{\mathrm{T}}{X}_f{)}^{-1}{X}_f^{\mathrm{T}}y\\ & +y^{\mathrm{T}}{X}_f(X_f^{\mathrm{T}}{X}_f{)}^{-1}{X}_f^{\mathrm{T}}{X}_f(X_f^{\mathrm{T}}{X}_f{)}^{-1}{X}_f^{\mathrm{T}}y\\ =& y^{\mathrm{T}}y-y^{\mathrm{T}}{X}_f(X_f^{\mathrm{T}}{X}_f{)}^{-1}{X}_f^{\mathrm{T}}y.\end{array}$$
And $$\begin{array}{rl}{X}_f^{\mathrm{T}}{X}_f& =\left[\begin{array}{ccc}1& \cdots & 1\\ \cos 2\pi f\cdot 1& \cdots & \cos 2\pi f\cdot n\\ \sin 2\pi f\cdot 1& \cdots & \sin 2\pi f\cdot n\end{array}\right]\left[\begin{array}{ccc}1& \cos 2\pi f\cdot 1& \sin 2\pi f\cdot 1\\ ⋮& ⋮& ⋮\\ 1& \cos 2\pi f\cdot n& \sin 2\pi f\cdot n\end{array}\right]\\ & =\left[\begin{array}{ccc}n& \sum _{t=1}^n\cos (2\pi ft)& \sum _{t=1}^n\sin (2\pi ft)\\ \ast & \sum _{t=1}^n{\cos }^2(2\pi ft)& \sum _{t=1}^n\cos (2\pi ft)\sin (2\pi ft)\\ \ast & \ast & \sum _{t=1}^n{\sin }^{2}2\pi ft\end{array}\right].\end{array}$$ Next let ${\displaystyle r={\mathrm{e}}^{2\pi \mathrm{i}ft}}$, ^[4]$$\begin{array}{rl}& \sum _{t=1}^n\cos (2\pi ft)=\sum _{t=1}^n\frac{{\mathrm{e}}^{2\pi \mathrm{i}ft}+{\mathrm{e}}^{-2\pi \mathrm{i}ft}}{2}\\ \stackrel{r={\mathrm{e}}^{2\pi \mathrm{i}f}}{=}& \frac{1}{2}\sum _{t=1}^n{r}^t+\frac{1}{2}\sum _{t=1}^n{r}^{-t}\\ =& \frac{1}{2}\frac{{\mathrm{e}}^{2\pi \mathrm{i}f}}{{\mathrm{e}}^{2\pi \mathrm{i}f}-1}({\mathrm{e}}^{2\pi \mathrm{i}f}-1)+\frac{1}{2}\frac{{\mathrm{e}}^{-2\pi \mathrm{i}f}}{{\mathrm{e}}^{-2\pi \mathrm{i}f}-1}({\mathrm{e}}^{-2\pi \mathrm{i}f}-1)=0.\end{array}$$

Next similarly $$\begin{array}{rl}& \sum _{t=1}^n{\cos }^2(2\pi ft)=\sum _{t=1}^n\frac{1+\cos (4\pi ft)}{2}=\frac{n}{2},\\ & \sum _{t=1}^n\cos (2\pi ft)\sin (2\pi ft)=\frac{1}{2}\sum _{t=1}^n\sin (4\pi ft)=0.\end{array}$$

So $$X_f^{\mathrm{T}}{X}_f=\left[\begin{array}{ccc}n& 0& 0\\ 0& \frac{n}{2}& 0\\ 0& 0& \frac{n}{2}\end{array}\right],(X_f^{\mathrm{T}}{X}_f{)}^{-1}=\left[\begin{array}{ccc}\frac{1}{n}& 0& 0\\ 0& \frac{2}{n}& 0\\ 0& 0& \frac{2}{n}\end{array}\right].$$ Then $$\begin{array}{rl}\mathrm{RSS}(f)& =y^{\mathrm{T}}y-(\sum _{t=1}^n{y}_t\sum _{t=1}^n{y}_t\cos 2\pi ft\sum _{t=1}^n{y}_t\sin 2\pi ft)\left[\begin{array}{ccc}\frac{1}{n}& 0& 0\\ 0& \frac{2}{n}& 0\\ 0& 0& \frac{2}{n}\end{array}\right]\left(\begin{array}{c}\sum _{t=1}^n{y}_t\\ \sum _{t=1}^n{y}_t\cos 2\pi ft\\ \sum _{t=1}^n{y}_t\sin 2\pi ft\end{array}\right)\\ & =y^{\mathrm{T}}y-\frac{1}{n}{\left(\sum _{t=1}^n{y}_t\right)}^2-\frac{2}{n}{(\sum _{t=1}^n{y}_t\cos 2\pi ft)}^2-\frac{2}{n}{(\sum _{t=1}^n{y}_t\sin 2\pi ft)}^2\\ & =\sum _{t=1}^n(y_t-\bar{y}{)}^2-\frac{2}{n}{(\sum _{t=1}^n{y}_t\cos 2\pi ft)}^2-\frac{2}{n}{(\sum _{t=1}^n{y}_t\sin 2\pi ft)}^2.\end{array}$$ Note again that this is only true when ${\displaystyle f\in [0,\frac{1}{2}]}$ and ${\displaystyle f}$ is Fourier frequency.

So $$\mathrm{RSS}(f)=\sum _{t=1}^n(y_t-\bar{y}{)}^2-2I(f).$$

We can also rewrite $$I(f)=\frac{1}{n}{|\sum _{t=1}^n{y}_t(\cos 2\pi ft+\mathrm{i}\sin 2\pi ft)|}^2=\frac{1}{n}{\left|\sum _{t=1}^n{y}_t{\mathrm{e}}^{-2\pi \mathrm{i}ft}\right|}^2.$$ This is exactly a Fourier transformation.

5 Some Other Nonlinear Regression Models

1. ${\displaystyle y_t={\beta }_0+{\beta }_{1}t+{\beta }_2\cos (2\pi ft)+{\beta }_3\sin (2\pi ft)+{\epsilon }_t}$. $$\mathrm{RSS}(f)=\sum _{t=1}^n(y_t-{\beta }_0-{\beta }_{1}t-{\beta }_2\cos (2\pi ft)-{\beta }_3\sin (2\pi ft){)}^2.$$
2. (Broken stick / change of slope) ${\displaystyle y_t={\beta }_0+{\beta }_{1}t+{\beta }_2(t-s{)}_{+}+{\epsilon }_t}$, here ${\displaystyle (t-s{)}_{+}=max\{t-s,0\}}$.